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A one species time-delay reaction-diffusion system defined on a complex networks is studied. 
Travelling waves are predicted to occur as follows a symmetry breaking instability of an homoge¬ 
nous stationary stable solution, subject to an external non homogenous perturbation. These are 
generalized Turing-like waves that materialize in a single species populations dynamics model, as the 
unexpected byproduct of the imposed delay in the diffusion part. Sufficient conditions for the onset 
of the instability are mathematically provided by performing a linear stability analysis adapted to 
time delayed differential equation. The method here developed exploits the properties of the Lam¬ 
bert W-function. The prediction of the theory are confirmed by direct numerical simulation carried 
out for a modified version of the classical Fisher model, defined on a Watts-Strogatz networks and 
with the inclusion of the delay. 
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Collective dynamics spontaneously emerge in a vast plethora of physical systems and often play a role of paramount 
importance for the efficient implementation of dedicated functions. Travelling waves are among the most studied 
phenomena for their ubiquitous and cross-disciplinary interest. Periodic traveling waves are for example encountered 
when describing self-oscillatory and excitable systems in different realms, from chemistry to biology, passing through 
physics. The Fisher 0,0 equation, introduced to characterize the spatial spread of an advantageous allele, defines 
the paradigmatic arena for addressing the peculiarities of travelling wave solutions in a reaction diffusion system 
in a spatial continuous domain. This is a one species population dynamic model, which assumes a logistic rule of 
replication and growth. The microscopic entities belonging to the scrutinized population can also delocalized in space, 
as follows a standard diffusion mechanism. The interplay between the aforementioned processes yields stable travelling 
wave solutions, which are selectively generated starting from a special class of initial conditions 0. More generally, 
it is however interesting to speculate on the possibility for a system to yield self-organized collective patterns of the 
travelling wave type, as follows a symmetry breaking instability, seeded by diffusion. Starting from a homogeneous 
solution subject to a tiny, non homogeneous, initial perturbation, a reaction diffusion system can destabilize via 
a dynamical instability, identified by Alan Turing in a seminal work 0. The Turing instability, as the process is 
nowadays called, can drive the emergence of non linear stationary stable patterns, if at least two species, the activator 
and inhibitors, are diffusiong and mutually interacting in the embedding environment. Alternatively, the Turing 
mechanims can instigate travelling wave solutions, provided at least three species, one of which mobile, are assumed 
to interact via apt non linear couplings. For a reaction diffusion system hosted on a discrete heterogeneous spatial 
support, namely a network, the instability can eventually set in for a two species model, if just the inhibitor is allowed 
to crawl from one node to its adjacent neighbors 0. Also in this case, three coupled species are the minimal request 
for a travelling wave to rise from a stochastic perturbation of an initial homogeneous stationary stable state. 

Starting from these premises, the aim of this Letter is to tackle the above problem under a radically different angle 
and thus introduce the simplest mathematical setting for which travelling wave solutions are generated, on a network, 
as a symmetry breaking instability of the Turing type. To anticipate our finding, we will prove that a one species 
model endowed with a delay in the diffusion can produce the sought instability. The analysis holds in general, but to 
demonstrate our conclusion we shall refer to a Fisher equation defined on a complex networks and modified with the 
inclusion of the delay term. 

The usage of time-delay differential equations (DDEs) defined on complex networks 0-0 is nowadays very popular, 
from pure to applied sciences [i®. Delays are for instance introduced to model finite communication or displace¬ 
ment time of quantities across network links. The imposed time-delay can non trivially interfere with the reactive 
dynamics, taking place on each node of the graph, thus resulting in unexpected emergent properties. For example, 
the classical paper m deals with time-delayed system made of two coupled phase oscillators and demonstrates the 
existence of multistability of synchronized solutions: an invariant manifold exists which attracts all the solutions of 
the system, yielding global oscillatory phenomena. Since this pioneering contribution, the subject has gained lot 
of attention and several contribution have been reported, assuming linear systems 121, coupled oscillators (l3l - [l5l| . 
oscillations death ( 161-00 and oscillations control [ 10 . Theory has been fruitfully applied to tackle real problems, see 



2 


for instance [20l - l^ . 

The work of this Letter moves from this reference context to build an ideal bridge with the Turing-like theory 
of pattern formation on complex networks. The discrete Laplacian operator, that encodes for the diffusion of the 
mobile populatiom incorporates a delay term. Our analysis of the DDEs exploits the properties of the Lambert 
W-function [^, thus the solution of the linearised equation can be given in closed analytical form. Furthermore, 
we have full access to the associated eigenvalues and their dependence on the involved parameters is made explicit. 
This allows us to go beyond the technique based on the computation of the Hopf bifurcation, namely to determine 
the parameters for which the eigenvalue with the largest real part passes through a pure imaginary value. 

We consider an undirected connected network composed by n nodes and assume one species to diffuse, from node 
to node, via the available links. Reactions also take place on each node, as dictated by a specific non-linear function 
/ of the local species concentration. Let us denote by Xi(t) the species concentration on node i, at time t. Then its 
time evolution is governed by 


ii{t) = - Tr)) + - Td) , (1) 
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where D stands for the diffusion coefficient, Lij = Gij — kiSij is the Laplacian matrix of the network whose adjacency 
matrix is given by G and ki = Gij identifies the degree of the i-th node, is the delay involved in the reaction 
occurring at each node, while Td is the delay due to the displacement across nodes. For the sake of simplicity we 
hereby assume the delay to be independent from the link indexes [3lj |. As mentioned earlier, it is well known that the 
above one species reaction-diffusion system cannot exhibit Turing-like instability, in the limiting case for Tr = Td = 0. 
As we shall argue, the introduction of a hnite delay = r > 0, will significantly alter this conclusion. 

To proceed in the analysis we assume a stable homogeneous equilibrium, Xi{t) = x for all i = 1,... ,n and t > 0 
and look for sufficient conditions to destabilize such equilibrium, as follows the introduction of a non homogeneous 
perturbation, which in turn activates the diffusion part. To determine the preliminary conditions that have to be 
met for the homogeneous equilibrium to be stable, we linearize system (ED with D = 0, around x, and recall that 
Tr = T > 0. Let A = f'{x), then the characteristic equation reads A = Ae~^'^ whose solutions are: 

Afe = —Wfc(rA), fc € Z, (2) 

T 

Wk being the k-th branch of the Lambert W-function To guarantee the needed stability of the homogeneous 

equilibrium x, one has to require that (A, r) and an integer k exist for which lRAfc(A, r) < 0. 

The Lambert W-function is the complex multivalued function of the complex variable z G C defined to be solution 
of the equation z = W(z)e^''^\ It has infinitely many branches denoted by Wk{z), for fc G Z; among them Wo{x) 
- the principal branch - is obtained by restricting z to lie on the real axis, more precisely on 5Rz G (—e“^, -boo), and 
with the constraint Wo(z) > — 1. Hence the branch cut of Wq is defined by {z : — oo < IRz < — , 3z = 0}. Let us 

observe that k = 0 and k = —1 are the only branches for which the Lambert W-function can assume real values. For 
all remaining fc, SITfc(z) ^ 0 for all z (see Fig. [T|). 

Roughly speaking, Wq bends the z plane (cut along IRz < —e“^) into a parabolic like domain in the w plane, whose 
boundary curves Sw i—>■ IRw = cotan Sw (blue solid and dotted curves in Fig. [T|) are bounded by tt and —tt. 
Moreover Wo satisfies the following relevant condition (Lemma 3 of [^ 1 

Vz G C : max5RWfc(z) = 5RWo(z). (3) 

The previous Eq. m allows to restate the stability condition of the homogeneous equilibrium as follows: 

3(A,t) such that 5RWo(rA) <0. (4) 

Hence, from the properties of the Lambert W-function represented in Fig. [T] this amounts to require 

-^<tA<{). (5) 

We now turn to considering the effect of a non homogeneous perturbation, superposed to the postulated homogenous 
equilibrium. This implies studying the full system (IB1|) . with Tr = Td = t > 0 and D > 0. Linearising as before, we 
get 


5xi{t) = A6xi{t — t) + D Lij6xj{t — t) . 

3 


( 6 ) 
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FIG. 1: The Lambert W function: principal branch Wo(z). Left panel: in the complex plane 2 G C, we represent the 
upper part of the branch cut {z : —00 < ?Siz < —e~^ , Qz = O'*"} by a solid line and the lower part of the branch cut 

{z : —00 < K 2 < —e~^ , Qz = 0~} by a dashed line; the circle denotes the point ( —l/e,0) while the square refers to (—7r/2,0). 

Right panel: the complex plane w € C, where w = Wo{z). The solid blue line is the image of the upper part of the branch 
cut Qw I—Kui = — Sju) cotanSrro for 0 < SJw < tt, while the dashed blue line is the image of the lower branch cut via Wo, 

Qw I—Km = —SrwcotanSrw for —tt < Qw < 0. The circle at coordinates ( — 1,0) is the image of the point (—l/e,0) and 

the squares positioned at (0,7r/2), respectively at (0, —7r/2), are the image of the point of coordinates (—7r/2,0^), respectively 
(—7r/2, 0~). The red dashed line is the image of the positive real axis while the green curved line is the image of the imaginary 
axis K 2 = 0. Observe that if = 0 and —n/2 < K 2 < —1/e, then —1 < ^Wo{z) < 0, |S>lFo(z)| < n/2 and S5Wo(2) / 0 
(blue solid and dotted lines), if $52 = 0 and — 1/e < K 2 < 0, then — 1 < ^Wo{z) < 0 and SJVFo(2) = 0 (blue dashed line), if 
^2 = 0 and K 2 < —7r/2, then ^Wo{z) > 0 and |£iVFo(2)| > 7r/2 (blue solid and dotted lines) and if Qz = 0 and K 2 > 0, then 
^Wo{z) > 0 and S>Wo(2) = 0 (red dashed line). 


The Laplacian matrix has a complete set of orthonormal eigenvectors associated to the topological eigenvalues 
0 = > ... > A", J2j for a = 1,... ,n. Using this basis to decompose Sxi{t) and employing 

once again the ansatz of exponential growth 

(5xi(t) = ^ for all (7) 

Oi 

we eventually get from (|6l) the characteristic equation 

Xo,-{A + DA°‘)e-^-^ = 0. ( 8 ) 

That is rAae'^'*'" = t{A + I1A“), whose solutions are: 

Xak = —Wk {t{A + DA°‘)), fc G Z and a = 1,..., n. (9) 

’ T 

Symmetry breaking instabilities seeded by diffusion can set in if at least a pair k and a exists such that > 0, 

provided — f < tA < 0 (the homogeneous equilibrium must be stable as a prerequisite of the analysis) [s^. When 
a bounded family of k exist for which KAa^fe > 0, the one with the largest real part dominates the instability and 
shapes the emerging pattern. Recalling again Eq. we finally get the following sufficient condition for the onset of 
Turing-like instability in a one species DDE of the general type dHU: 

da G [2,..., n] such that —KWq {t{A + DA°‘)) > 0 . (10) 

T 

Exploiting the properties of the Lambert W-function represented in Fig. [U Eq. m is satisfied whenever [sj there 
exists d > 1 such that 

TDA“<-|-rA. (11) 

Observe also that for such a, 3Ao,a 7 ^ 0. Hence, the instability materialize in the appearance of travelling waves. 
Since KWo(a;) is increasing, for decreasing real x < —tt 12, the dispersion relation attains its maximum at a = n, the 
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Laplacian eigenvalues being ordered for decreasing real parts. Hence, the Turing-like instabilities cannot develop if 
tDK^ > Summing up, travelling waves are expected to develop, for a fixed network topology, and sufficently 

large values of r and D. Conversely, for a fixed choice of the parameters (t, D and H), one should make the network 
big and so force A” to be large enough, in absolute value. 

As our goal is to determine the minimal model for which the aforementioned instability sets in, we now consider 
the reduced case for = 0 and = r > 0, that is the delay is only associated to the diffusion part. For D = Q the 
homogeneous equilibrium i, is stable if and only if A < 0. By repeating the procedure highlithed above we end up 
with the following characteristic equation: 

Aa - A - i:iA“e-^“^ = 0, (12) 

whose solution reads 

= ivFfe(riAA“e-"^) + A, fc G Z. (13) 

r 

The equilibrium x is hence destabilised by diffusion if there exist a > 1 for which 

^Wo{TDA°‘e-^^) > -At. (14) 

Let us observe that —At is positive and TDA°‘e~'^^ is negative. Moreover, we already remarked that ^Wo{x) is 
positive and increasing for decreasing negative x. Hence, a critical Xc < 0 exists [s^ for which lR.Wo{xc) = —At and 
Eq. (fTTl) is satisfied for all TDK°‘e~'^^ < x^- Travelling waves are hence predicted to manifest, for sufficiently large 
A“, in absolute value. We again stress that 9Ao,a ^ 0. 

To complete the general discussion we consider the dual problem, where Tr = t > 0 and t^ = 0. As already 
observed, the homogenous equilibrium x is stable if —7r/2 < rA < 0. The characteristic equation associated to this 
problem can be cast in the form: 

A« - Ae-^“^ - Z)A“ = 0, (15) 


whose solution is 

A„fc =iAA“ + iWfc(TAe-"^^“), fcGZ. (16) 

T 

One can prove that DA°‘ + ^Wo{TAe~^^^°‘) < 0 for all A“, a > 1. We are consequently led to conclude that the 
stable homogeneous equilibrium cannot undergo a diffusion driven Turing-like instability, if the delay term is solely 
confined in the reaction part. 

As an application of the previous theory, we take / in Eq. dHU to be logistic function /(x) = ax(l — x), as in 
the spirit of the Fisher model. At variance with the Fisher equation P, 0, we now imagine the species to be hosted 
on a discrete support, rather than a continuum segment. In the original Fisher scheme the emerging wave relates 
to the heteroclinic orbit of the system and requires a specific, step-like, initial profile. In our case the travelling 
wave will originate as follows a symmetry breaking instability of an initial randomic perturbation. To carry out the 
analysis, we select the homogeneous equilibrium solution x = 1 and look after to its associated stability properties, 
as a function of a, r, D and the network topology. Silencing the diffusion, ZI = 0, the stability condition Eq. 
rewrites 0 < ar < ^. In demonstrating our findings, we consider a Watts-Strogatz network made of 100 nodes, 
with average degree (fc) = 6 and probability to rewire a link p = 0.03. Other network topologies can be in principle 
assumed, returning similar qualitative conclusions. The chosen network is large enough so to have one eigenvalue for 
which the dispersion relation has positive real part (see left panel of Fig. [3]): thus Turing-like waves do exist. In Fig. [2] 
we report the result of a numerical solution of the system under scrutiny obtained with a RK4 method adapted to 
deal with delay differential equation. The parameters are set to the values a = 1.4, D = 0.05 and r = 1. The DDE 
should be complemented with the value of the function on the delay interval [—t, 0). In the spirit of a perturbation 
of the stable equilibrium, we decided to set Xi(t) = \-\- 5i for all t G [—r, 0) where 5i are random Gaussian numbers 
drawn from Af(0, 0.01). Observe that ar = 1.4 < -k/2 and thus the equilibrium x = 1 is stable in absence of diffusion. 
On the other hand, one can clearly appreciate that after a transient period, patterns do manifest as stable oscillations 
around the solution x = 1. 

In Fig. [3] the dispersion relation are displayed. In the left panel the quantity 

maxSJAa k = o = — 3?Wo (—ra -I- tDA°‘) , 

k ’ ’ T 

is plotted as a function of A“, for fixed a, D and r (as specified in Fig. [5]), and for the same Watts-Strogatz network. 
On can clearly identify several eigenvalues, for which 5 RAq,^o > 0- The right panel reports the maximum of the 
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FIG. 2: Numerical solution of the system under scrutiny with parameters a = 1.4, D = 0.05 and r = 1. Initial conditions are 


set equal to Xi{i) = 1 + (5i for all t £ 
network is a Watts-Strogatz network 
p = 0.03. 


—r, 0) where 5i are random gaussian numbers drawn from ^"(0, 0.01). The underlying 
l30l| made of 100 nodes, with average degree (fc) = 6 and probability to rewire a link 


dispersion relation as a function of the eigenvalues, i.e. maxa 5RAct,fc, versus (a, r) for fixed D = 0.05 and for the 
same Wattz-Strogatz network. The stability domain of the equilibrium without diffusion is bounded by a > 0 (r > 0, 
for physical reasons) and ar < 7r/2 (dash-dotted black curve). The Turing-like waves can emerge for all pairs (a,r), 
inside such limited domain, for which maxo, JRA^.o > 0. 



A“ 



FIG. 3: Dispersion relation. Left panel KAc.o as a function of A“, for fixed a, D and r, as specified in Fig. [2] and for the 
same Wattz-Strogatz network. Blue dots represent the values of A^.o computed for a given A“, while the solid black line is the 
continuous approximation. The horizontal dotted black line stands for the 0-th level. Right panel: max^ KA^.o, as a function 
of (a, r), for fixed D = 0.05 using the same Wattz-Strogatz network as employed in Fig. [2] The dash dotted curve ar = 7r/2 
delineates the boundary (together with a > 0 and r > 0) of the stability region of the homogenous solution x = 1. Hence, all 
pairs (a,r) for which 0 < ar < 7r/2 correspond to stable solution of the homogeneous equilibrium. Turing-like instabilities are 
thus allowed to develop for all pairs (a, r) inside such domain, for which max^ KAa,o > 0. The star refers to the setting of the 
simulation reported in Fig. [2] 

In conclusion, we have here shown that a one species time-delay reaction-diffusion system defined on a complex 
networks can exhibit travelling waves, as follows a symmetry breaking instability of an homogenous stationary stable 
solution, subject to an external non homogenous perturbation. These are Turing-like waves which emerge in a minimal 
model of single species population dynamics, as the unintuitive byproduct of the imposed delay. Based on a linear 
stability analysis adapted to time delayed differential equation, we provided sufficient conditions for the onset of 
the instability, as a function key quantities, as the reaction parameters, the delay, the diffusion coefficient and the 
network topology. The wave possesses multiple fronts and persists in time, without fading away as it happens for 
the customary Fisher equation. The observation that Turing-like instability can originate for a one species model 
evolving on a heterogeneous graph, provided a delay is included in the transport term, enables us to significantly 
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relax the classical constraints for the patterns to emerge and opens up the perspective for intriguing developments in 
a direction so far unexplored. 

The work of J.P, T.C. and M.A. presents research results of the Belgian Network DYSCO (Dynamical Systems, 
Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian 
State, Science Policy Office. 


Appendix A: Turing-like instability when the delay is confined in the diffusion term. 

Let us consider the case where the delay is present only in the diffusion part = t > 0 and = 0) 

i*(t) = f{x^{t))+D^L,jXj{t-T). (Al) 


To study the stability of the equilibrium x in this case we have to look for a > 1 for which 

-At. (A2) 

Let z = x + iy and w = ^ + iri, then the following relation is recovered once we impose z = we’", namely w = Wo( 2 ): 


X = cos rj — rj sin rf) 
y = e^{r] cos y + ^ sin 77 ). 


(A3) 


Let us introduce Xc = tDA°‘ and = —tA, values for which the equality holds in (IA2I) . Observe that Xc < 0 and 
> 0. Rewriting Eq. (IA3I1 for Zc = Xc (namely j/c = 0) and Wc = ^c + iVc we get: 

\xc = 6^“= (^c COS rjc - 7]c sin rjc) 

10 = ( 77 c cos 77c + Cc sin 77 ^). 

From the second equation one can obtain (implicitly) rjc as a function of ^c- 


rjc = -ic tan rjc ■ 


Because < 0 one can find a unique solution rjci^c) G ('^/2, t) (and the opposite one). Inserting this result into the 
first equation we get 


^c(^c) 


g-» 7 c(Cc)cotan ?7c(5c) 


Sm77c(^c) ’ 


So in conclusion given —tA one can obtain the critical Xc(—tA) for which we have equality in Eq. (IA2p . and thus 
conclude, by invoking the scaling properties of Wqi that for all A“ < Xc{—tA)/{tD) the strict inequality sign holds 
in Eq. 


Appendix B: Turing-like instability are impeded when the delay only appears in the reaction term. 

Let us consider system 


Xiit) = f{xi{t - t)) + D^L^jXj{t ), (Bl) 

3 

once the delay dependence is given only through the reaction term, Ti- = t > 0 and 
stability of the homogeneous equilibrium without diffusion is given by — ^ < tA < 0 . 
outlined in the main body of the paper - linearising, then expanding the perturbation 
of the Laplacian - one obtains the following characteristic equation: 

Aa - - DA“ = 0, (B2) 

that can be solved using the W-Lambert function to give 

= DA“ + -WkirAe-^^^'^), k€Z. (B3) 

T 


Td = 0. The condition for the 
Following the same procedure 
in the basis of the eigenvectors 
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Turing-like instability can emerge if the homogeneous equilibrium becomes unstable in presence of the diffusion. 
We then look for a and k such that SRAq ^ > 0. Using Lemma 3 of [2^ this is equivalent to 

> -TDh °‘, 

Let us observe that —tDA‘^ > 0. Because tA < 0, one cannot solve the previous equation with 'AWo{TAe~'^^^°‘) = 0 
(in this case one should have the argument of Wq to be positive). Hence we look for TAe~'^^^ < —ti 12. 

Let us introduce s = —DK'^ > 0, u = tA G (— 7 r/ 2 , 0 ) and the function 

g(s) = —s -I- , (B4) 

our goal is to prove that g(s) < 0 for all s > 0 which is in turn equivalent to stating that Turing-like instability cannot 
develop. 

Let us rewrite Eq. (IA3I) for x = ue^ and y = 0: 

{ ite® = cos r] — rj sin 77 ) 

0 = (77 cos rj + ^ sin 77 ). 

Isolating ^ in the second equation and inserting it in the first one, we get: 


M sin 77 


One can thus rewrite g{s) as follows: 


g{s) = -s + WVoiue’') = -s -|- ^ = -.J - log ( -?—^ -I- ^ 

\ M sm 77 / 

= - log ( -. 

\ usmriJ 

Because u > — 7 r /2 and tt/2 < 77 < tt we obtain — 77 /( 77 sin 77 ) > 1 and thus log > 0. 
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